Introduction

Depression is one of the most common mental health conditions worldwide. It can have significant consequences for individuals, communities, and wider society, including inmpacts on wellbeing and employment. One form of depression known to worsen during the winter months is Seasonal Affective Disorder (SAD) (springer/references).

The Covid-19 pandemic created major distruptions to daily life, including prolonged periods spent indoors due to lockdowns. Reduced exposure to natural light is a known risk factor for SAD(NHS), and the combination of winter daylight shortages and enforced indoor isolation may have contributed to increased levels of depressive symptoms during the pandemic. These changes may be reflected in pattern of antidepressant prescribing, particularly during winter periods.

The aim of this study is to investigate wether the COVID-19 pandemic influenced antidepressant prescribing and whether there are inequalities or geographic differences associated with this trend.

Methods

All data used in this study was obtained from publicly available NHS Scotland and Scottish Government sources.Antidepressant prescribing data were downloaded from the Data by Prescriber Location dataset, focusing specifically on winter periods (December to Februray) for the years 2019-2023, from the Public Health Scotland OpenData portal. Winter population estimates were taken from the National Records of Scotland census release.Socioeconomic information was sourced from the Scottish Index of Multiple Deprivation (SIMD 2020v2) datasets published by the Scottish Government. Spatial analysis was supported by incorporating NHS Scotland Health Board boundary shapefiles, which are publicly available through PHS and used to map regional variation in prescribing.

Antidepressant medications were identified using the BNF item code prefix “0403”, and all prescriptions with codes starting with “0403” were retained for analysis. These filtered data were then processed and summarised to generate visualisations of prescribing patterns.

For the purpose of this study, the COVID-19 period is defined as March 2020 to June 2021, corresponding to the main period of national lockdowns. This places Winter 2019 as pre-COVID, Winter 2020/2021 as during COVID, and Winter 2022 as post-COVID.

##Loading Necessary Packages
library(tidyverse) 
library(here) # directory stucture
library(gt) # tables
library(janitor) # cleaning data
library(ggplot2) # plotting graph
library(sf) # to read in map data 
library(readxl) # to read in map data
library(plotly) # to make interactive
library(viridis) # colour palette for map
library(ggiraph)
# Loading multiple data files (18 CSVs) containing counts of medications dispensed in the community
files <- list.files(here("data", "winter_data"), pattern = "csv")
winter_data <- files %>% # Read all downloaded winter prescribing files stored in the 'winter_data' folder
  map_dfr(~read_csv(here("data", "winter_data", .))) %>% 
clean_names()#combine into one dataframe and clean column names
filtered_winter_data <- winter_data %>% 
filter(str_starts(bnf_item_code,"0403")) %>%  #Filter for antidepressants only (BNF item code starts with "0403" 
  mutate(year = as.numeric(substr(paid_date_month,1,4)), month = as.numeric(substr(paid_date_month,5,6))) %>% #Extract year and month from 'paid_date_month' for grouping into winter periods
  mutate(winter_year=case_when(month == 12 ~ year + 1, 
month %in% c(1,2) ~ year) )  #Create 'winter_year' column: December counts belong to the next year, Jan-Feb to current year
filtered_winter_data <- filtered_winter_data %>% 
  unite("healthboards",hbt2014,hbt,sep = "_") #Merge Health Board columns into a single column since they contain the same information
filtered_winter_data$healthboards <- gsub("[NA]","",filtered_winter_data$healthboards) 
    filtered_winter_data$healthboards <-
      gsub("_","",filtered_winter_data$healthboards)#Remove any remaining 'NA' strings and underscores to standardize Health Board codes

Figure 1

 #summarise total antidepressant prescriptions for each winter year
winter_years_data <- filtered_winter_data %>% 
  group_by(winter_year) %>% 
  summarise(total_items=sum(number_of_paid_items,na.rm = TRUE))
#plot line graph to visualise overall trend in total prescriptions across winter years
plot <- suppressWarnings( ggplot(winter_years_data,aes(x=winter_year,y=total_items)) + #line connecting total prescriptions per year
  geom_line(linewidth=0.7,colour = "blue") +
  geom_point(aes(text = paste0(
      "<b>Winter Year:</b> ",winter_year, "<br>",
      "<b>Total Items:</b> ",scales::comma(total_items)),
size = 1))+ #points with interactive tooltip showing year and total
  geom_vline(xintercept = 2020, 
linetype = "solid", colour = "brown") + geom_vline(xintercept = 2022, 
linetype = "solid", colour = "brown") + #reference lines to mark COVID-related winters
  scale_x_continuous(breaks=2017:2023) + #ensures all years appear on x-axis
  labs(title="Total Antidepressant Prescriptions During Winter Season",x="Winter Year",y="Total Antidepressant Prescriptions") +
  theme_minimal() )#thiscode suppresses warnings
  
ggplotly(plot,tooltip="text") #convert ggplot to interactive plotly object with tooltips

The line graph shows a steady upward trend in winter antidepressant prescribing across all years. Although the COVID-19 period (2020–2022) is highlighted, prescribing patterns do not markedly diverge from pre-pandemic trends, with increases during COVID appearing slightly less steep than earlier years. The y-axis is truncated to improve visibility of year-to-year changes. 

#Prepare Population
population <- readxl::read_excel(here("data","population.xlsx"), skip=10) %>% 
  clean_names() %>% 
  select(x2,all_people) %>% 
  filter(!is.na(all_people)) %>% 
  rename(h_bname = "x2",hb_population = "all_people") %>% 
filter(!str_detect(hb_population,"Cells"))

filtered3_winter_data <- filtered_winter_data %>% 
  group_by(healthboards,winter_year,gp_practice) %>% 
  summarise(paid_quantity = sum(paid_quantity,na.rm=TRUE)) 

SIMD <- readxl::read_excel(here("data","SIMD.xlsx")) %>% 
  clean_names() 

combined_pop_simd <- SIMD %>% 
  full_join(filtered3_winter_data,join_by(h_bcode == healthboards)) %>% 
  left_join(population,join_by(h_bname == h_bname))

Figure 2

#sum antidepressant prescriptions per health board and winter year (from 2019 onwards)
hb_prescribing <- filtered_winter_data %>%
  filter(winter_year >= 2019) %>%
  group_by(healthboards, winter_year) %>%
  summarise(total_paid = sum(paid_quantity, na.rm = TRUE), .groups="drop")
#clean SIMD data (keep unique health board codes and names)
simd_clean <- SIMD %>% select(h_bcode, h_bname) %>% distinct()
#combine with prescribing data with SIMD and population data,calculate prescriptions per head
combined <- hb_prescribing %>%
  left_join(simd_clean, join_by(healthboards == h_bcode)) %>%
  left_join(population, by = "h_bname") %>%
  mutate(quantity_per_head = total_paid / hb_population) %>%
  filter(!is.na(h_bname))
#read and clean NHS health board shapefile
NHS_healthboards <- st_read(here("data","Week6_NHS_HealthBoards_2019.shp")) %>%
  clean_names() %>% rename(h_bname = hb_name)
## Reading layer `Week6_NHS_healthboards_2019' from data source 
##   `/Users/olufimihanfaturoti/Year 3 Medicine/data-science/B251495/data/Week6_NHS_healthboards_2019.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 14 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 7564.996 ymin: 530635.8 xmax: 468754.8 ymax: 1218625
## Projected CRS: OSGB36 / British National Grid
#join spatial data with prescribing data, fill missing value with 0
retry_mapped_data <- NHS_healthboards %>%
  left_join(combined, by = "h_bname") %>%
  st_as_sf() %>%
  mutate(quantity_per_head = replace_na(quantity_per_head, 0))
#plot interactive map of antidepressant prescriptions per head
plot_map <- retry_mapped_data %>% 
  ggplot() +
  geom_sf_interactive(aes(fill = quantity_per_head,
        tooltip = paste0(h_bname,
                         "\nWinter Year: ", winter_year,
                         "\nItems per Head: ", round(quantity_per_head,2))),
    colour = "white",
    size = 0.1) +
  scale_fill_viridis_c(option = "mako", direction = -1, name = "Items per Head") + #blueish colour palette
  labs(title="Antidepressant Prescriptions per Head",
       subtitle="By Health Board and Winter Year") +
  facet_wrap(~winter_year) + # separate maps per year
  theme_void() +
  theme(strip.text = element_text(size=12, face="bold"),
        plot.title = element_text(face="bold", size=16),
        plot.subtitle = element_text(size=10))
#convert to interactive map
interactive_map <- girafe(ggobj = plot_map);interactive_map

The map shows the geographic distribution of antidepressant prescriptions per head, revealing a consistent north–south gradient with higher rates in central urban regions. Post-2020 intensification suggests a system-wide rise in prescribing, with areas such as Lanarkshire and Glasgow showing the largest increases compared with 2019. Regions with historically higher prescribing appear to have experienced the greatest growth during and after the pandemic.

Figure 3

#| fig-align: center
#| fig-width: 6
#| fig-height: 6
#| fig.cap: "The boxplot shows no strong socioeconomic gradient in antidepressant prescribing across SIMD quintiles, with median rates remaining relatively flat in both 2019 and 2022. Although overall prescribing increased uniformly after COVID-19, greater variability and more extreme outliers in 2022 suggest widening differences between individual health boards or practices, rather than between deprivation groups."
#| message: false


winter_population <- list.files(here("data", "winter_population"), pattern = "csv", full.names = TRUE) %>%
  map_dfr(read_csv) %>%
  clean_names() %>%
  select(date, practice_code, hb, all_ages) %>%
  mutate(winter_year = as.numeric(substr(date, 1, 4)),
    winter_year = case_when(
      winter_year == 2020 ~ 2019,
      winter_year == 2023 ~ 2022,
      TRUE ~ winter_year ),
    healthboards = hb) %>%
  select(-hb)
#summarise and filter combined population data for selected years
summary_winter_with_simd <- combined_pop_simd %>%
  filter(!is.na(simd2020v2_quintile), winter_year %in% c("2019","2022")) %>%
  group_by(simd2020v2_quintile, h_bname, gp_practice, winter_year) %>%
  summarise(avg_paid_over_winters = mean(paid_quantity), .groups = "drop") %>%
  rename(practice_code = gp_practice) %>%
  mutate(winter_year = as.numeric(winter_year))
# Join datasets
joined_data <- summary_winter_with_simd %>%
  full_join(winter_population, by = c("winter_year", "practice_code"))
# Calculate prescriptions per 1000
final_boxplot <- joined_data %>%
  group_by(simd2020v2_quintile, h_bname, practice_code, winter_year) %>%
  reframe(
    avg_per_1000 = (avg_paid_over_winters / all_ages) * 1000,
    year_label = factor(winter_year,
                        levels = c(2019, 2022),
                        labels = c("Pre-COVID (2019)", "Post-COVID (2022)"))
  )

#create boxplot of average prescriptions per 1000 by SIMD quintile, faceted by pre/post COVID
box2 <- final_boxplot %>%
  ggplot(aes(x = simd2020v2_quintile, y = avg_per_1000, fill = factor(simd2020v2_quintile))) +
  geom_boxplot() +
  scale_fill_viridis(discrete = TRUE, alpha = 0.6) +
  scale_y_continuous(labels = scales::label_comma()) +
  coord_cartesian(ylim = c(0, 70000))  + #scaled y axis graph to be able to visualise better
  theme_minimal() +
  theme(legend.position = "none") +
  facet_wrap(~year_label) +
  labs(
    title = "Average Winter Antidepressant Prescriptions Per by SIMD Quintile",
    x = "SIMD Quintile (1 = Most Deprived)",
    y = "Avg prescriptions per 1000"
  )

ggplotly(box2)

Figure 4

#| fig-align: center
#| fig-width: 6
#| fig-height: 6
#| fig.cap: "The lollipop chart highlights strong geographical inequalities in the impact of COVID-19 on antidepressant prescribing. Urban health boards(Lanarkshire and Greater Glasgow and clyde) show increases up to ten times greater than rural and island areas, suggesting that factors such as population density, service pressure, and differing pandemic experiences played a larger role than deprivation alone. message: false
#create population-weighted deprivation profile per health board
hb_deprivation <- SIMD %>%
  group_by(h_bname) %>%
  summarise(#weighted mean SIMD quintile (lower value = more deprived)
    avg_quintile = weighted.mean(simd2020v2_quintile, w = population, na.rm = TRUE),
#percentage of population living in most deprived quintiles (Q1 & Q2)
    pct_deprived = sum(population[simd2020v2_quintile <= 2], na.rm = TRUE) /sum(population, na.rm = TRUE) * 100,
#total population per health board (for reference)
  total_pop = sum(population, na.rm = TRUE)
  ) %>%
  arrange(avg_quintile)  #rank health boards from most → least deprived


#join deprivation ranking to prescribing data (2019 = pre-COVID, 2022 = post-COVID)
final_lollipop_data <- combined %>% 
  filter(winter_year %in% c("2019", "2022")) %>%      #keep only pre and post pandemic years
  mutate(period = ifelse(winter_year == 2019, "pre", "post")) %>%  #label periods
  left_join(hb_deprivation, by = "h_bname")            #add deprivation profile
#calculate % change in prescribing per health board
lollipop_plot <- final_lollipop_data %>% 
  group_by(h_bname, period) %>% 
  summarise(
    avg_prescribe = mean(total_paid),     
    avg_quintile = first(avg_quintile), 
    .groups = "drop"
  ) %>% 
  pivot_wider(
    names_from = period, 
    values_from = avg_prescribe           #create 'pre' and 'post' columns
  ) %>%
  mutate(
    pct_change = round(((post - pre) / pre) * 100, 2)  # Percentage change
  ) %>%
  arrange(avg_quintile) %>%               #rank health boards by deprivation
  mutate(h_bname = factor(h_bname, levels = unique(h_bname)))  
#lollipop Plot (with hover labels showing HB name + % change)
pct_lollipop <- lollipop_plot %>% 
  ggplot() +
  geom_segment(
    aes(
      x = h_bname, xend = h_bname,
      y = 0, yend = pct_change,
      # Hover tooltip content for stems
      text = paste0(
        "Health Board: ", h_bname, "<br>",
        "Change: ", pct_change, "%" )),
    color = "skyblue" ) +
  geom_point(
    aes( x = h_bname, y = pct_change,
      text = paste0(
        "Health Board: ", h_bname, "<br>",
        "Change: ", pct_change, "%")),
    color = "blue", size = 3, alpha = 0.6)+
  theme_light() +
  coord_flip() +
  theme_minimal()+
  labs(
    title = "Percentage Change in Antidepressant Prescribing per Health Board (Pre vs Post COVID-19)",
    x = "Health Board (ranked from most deprived to least deprived)",
    y = "Percentage Change (%)"
  )

# Convert to interactive plot with hover labels
ggplotly(pct_lollipop, tooltip = "text")

The lollipop data provides a summary of the impact of covid 19 . overall each health board experienced and increase with deprevation not having as much of an impact as expects.However lankashire and glasgow has the bggest increase which could potential lbe linked to be more deprived or larger population.

Conclusion Limitations Next Steps